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ABSTRACT 


During  this  contract’s  performance  period  (1  October  2004-30  September  2007)  we  intend  to  complete  the 
development  of  a  method  for  determining  crust  and  upper  mantle  structure  in  aseismic  regions  and  evaluate  its 
capacity  to  produce  sufficiently  accurate  models  for  locating  hypocenters  of  small-to-moderate  seismic  events. 
Specifically,  we  will 

•  complete  development  of  a  waveform  modeling  approach  to  determine  the  structure  of  the  crust  and  upper 
mantle  in  a  broad  region  surrounding  a  broadband  seismographic  station; 

•  validate  the  method  through  careful  comparison  between  models  produced  with  our  method  and  models 
produced  by  higher  resolution  active-source  experiments,  as  well  as  models  produced  by  surface  wave 
inversion  and  receiver  function  methods; 

•  apply  and  validate  methods  for  computing  uncertainty  statistics  from  the  products  of  global  optimization 
searches;  and 

•  refine,  apply,  and  evaluate  waveform  correlation  methods  to  ground-truth  data  from  regional  events  and 
explosions  with  well-constrained  focal  depths,  using  crust/upper  mantle  models  produced  with  our  method; 

•  apply  the  modeling  method  to  teleseismic  data  recorded  at  stations  in  China,  Africa,  and  the  Middle  East; 

•  apply  and  evaluate  waveform  correlation  methods  to  determine  focal  depths  for  small-  and  moderate- 
magnitude  regional  events  recorded  at  stations  in  China,  Africa,  and  the  Middle  East  using  crust/upper 
mantle  models  determined  in  this  study. 

The  high  frequencies  and  long  time-series  required  for  phases  that  arrive  near  the  direct  SV  phase,  including  Sp, 
SsPmP,  and  shear-coupled  PL  waves,  their  deep  penetration  into  the  Earth  and  observation  at  teleseismic  distances 
make  the  computation  of  synthetic  seismograms  time-consuming.  We  have  parallelized  and  optimized  a  synthetic 
seismogram  code  based  on  the  reflectivity  method  and  are  now  able  to  compute  complete  seismograms  up  to  0.5  Hz 
in  just  over  one  minute  using  eight  AMD  Opteron  processors.  This  has  allowed  us  to  produce  the  best- fitting 
waveform  matches  that  result  from  searches  of  many  thousands  of  models,  yet  searching  the  model  space 
comprehensively  multiple  times,  as  is  required  to  estimate  posterior  covariance,  correlation,  and  posterior 
probability  distribution  functions  is  still  prohibitive.  However,  the  speed-up  in  computation  time  is  nearly  linear 
with  the  number  of  processors  used,  so  a  massively  parallel  computer  installed  recently  at  the  University  of  Texas  at 
Austin  should  allow  us  to  compute  posterior  probability  density  function  (PPD)  estimates  readily.  We  intend  to 
explore  the  utility  of  PPD  estimates  for  estimating  the  reliability  of  our  crustal  models  for  the  purposes  of  estimating 
hypocenters  of  small-  to  moderate-sized  regional  seismic  events. 
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OBJECTIVES 

Our  research  focuses  on  developing  and  testing  methods  for  determining  event  locations  at  regional  distances  (100- 
1 ,000  km)  using  one  or  a  few  three-component,  broadband  seismographic  stations.  To  be  useful  for  discriminating 
explosions  from  earthquakes,  focal  depths  must  be  determined  accurately  to  within  a  few  kilometers,  yet  the 
structure  of  the  crust  and  upper  mantle  between  the  source  and  station  strongly  influences  the  arrival  times  and 
amplitudes  of  regional  phases.  Our  first  task,  therefore,  is  to  determine  this  structure.  We  will  attempt  to  obtain  a 
sufficiently  accurate  model  of  crustal  structure  beneath  a  given  station  by  searching  over  a  wide  variety  of  crustal 
models  to  find  the  one  whose  synthetics  best  match  the  amplitudes  and  travel  times  of  phases  that  arrive  in  the  time 
window  around  the  direct  S  phase,  and  which  appear  most  strongly  on  the  radial  component  seismogram.  We  use 
data  from  large-magnitude  (6<Mb<7)  deep  events  located  at  teleseismic  distances  for  this  purpose,  for  reasons 
discussed  in  detail  below. 

After  structural  models  are  in  hand  we  propose  to  find  focal  depths  by  means  of  a  waveform  correlation  technique. 
We  will  compute  suites  of  synthetic  seismograms  for  the  candidate  model  corresponding  to  a  given  broadband 
station  and,  given  data  from  small-magnitude  regional  events  recorded  at  the  same  station,  compute  correlations 
between  the  suites  of  synthetic  seismograms  generated  for  a  variety  of  distances  and  focal  depths. 

Case  studies  around  the  world,  including  application  to  “ground-truth”  data  in  which  focal  depths  are  well- 
constrained,  are  required  in  order  to  determine  the  range  of  applicability  and  usefulness  of  our  event  location 
strategy.  Small  magnitude  events  (2<mb<4)  are  rarely  recorded  at  more  than  a  few  stations  in  many  parts  of  the 
world,  yet  identifying  these  events  is  essential  to  monitoring  nuclear  explosions  effectively.  Our  objective,  therefore, 
is  to  develop  and  test  methods  for  distinguishing  between  explosions  and  natural  events  in  sparsely-instrumented 
parts  of  the  world.  Our  strategy  is  to  estimate  hypocenter  locations  for  small-  and  moderate-magnitude  events  using 
one  or  a  few  single  broadband,  three-component  stations.  Focal  depth  is  a  particularly  helpful  discriminating 
characteristic,  if  determined  reliably. 

Historical  Perspective 

The  problem  of  finding  earthquake  locations  using  only  a  single  three-component  station,  or  a  very  sparse  network 
of  stations,  has  received  a  great  deal  of  attention.  Some  authors,  recognizing  the  greater  difficulty  in  constraining 
focal  depths,  have  focused  on  estimating  event  epicenters  and  origin  times  (e.g.,  Magotra  et  al.,  1987;  Roberts  et  al., 
1989;  Kedrov  and  Ovtchinnikov,  1990;  Kim  and  Wu,  1997).  Others  have  attempted  to  estimate  not  only  focal 
depths  but  focal  mechanisms,  as  well,  via  waveform  modeling  of  regional  events  (Jimenez  et  al.,  1989;  Fan  and 
Wallace,  1991;  Zhao  and  Helmberger,  1994;  Walter,  1993;  Zhu  and  Helmberger,  1996).  Frohlich  and  Pulliam 
(1999)  review  efforts  to  locate  earthquakes  using  a  single  three-component  station  in  detail.  In  contrast,  the  problem 
of  determining  earthquake  focal  depths  with  one  or  a  few  stations  has  received  far  less  attention,  yet  obtaining 
accurate  estimates  of  focal  depths  is  important  to  tectonic  interpretations  of  seismicity,  to  understanding  seismic 
hazard,  and  to  seismic  monitoring  of  underground  nuclear  tests  (National  Research  Council,  1997).  If  a  seismic 
event  were  known  reliably  to  have  occurred  at  a  depth  greater  than  a  few  kilometers,  one  could  confidently 
categorize  that  event  as  an  earthquake  rather  than  an  explosion. 

The  most  common  approaches  to  determining  focal  depths  utilize  travel  times  (e.g.,  Douglas,  1967),  which  must  be 
picked  for  the  major  seismic  phases  and  then  back-projected  by  a  nonlinear  or  bootstrapped  linear  algorithm.  A 
great  deal  of  effort  has  been  devoted  to  methods  for  picking  arrival  times  automatically  (Roberts  et  al.,  1989;  Saari, 
1991 ).  While  these  methods  have  proven  useful  for  moderate  to  large  magnitude  events  at  far-regional  and 
teleseismic  distances,  particularly  for  events  that  are  deeper  than  a  few  tens  of  kilometers,  they  are  prone  to  picking 
errors.  These  picking  errors  become  relatively  more  important  and  more  problematic  for  location  procedures  when 
dealing  with  small  magnitude  and  shallow  events.  Furthermore,  shallow  events  have  little  time  separation  between 
the  downgoing,  direct  body  wave  phases  and  the  upgoing,  reflected  (“depth”)  phases  that  are  most  useful  for 
constraining  focal  depth.  One  must  use  high-frequency  data  to  have  any  hope  of  identifying  and  picking  distinct 
arrivals  for  these  phases,  which  again  complicates  the  picking  process  and  increases  the  likelihood  that  errors  will 
contaminate  the  location  process.  In  some  parts  of  the  world,  small-magnitude  and  relatively  shallow  events 
observed  at  regional  distances  typically  have  emergent  rather  than  impulsive  first  arrivals,  which  renders  travel  time 
picking  even  more  prone  to  errors.  Lastly,  reliable  estimates  of  both  epicenters  and  focal  depths  using  travel  times 
require  redundancy,  i.e.,  at  least  several  and  preferably  many  recordings  from  stations  that  are  well-distributed  with 
respect  to  azimuth  around  the  event.  But  the  great  majority  of  earthquakes  are  small-magnitude,  shallow  events, 
which  are  much  more  likely  than  large  events  to  be  recorded  by  just  one  or  a  few  seismographic  stations.  Locating 
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these  small  events  accurately  can  be  most  useful  for  discriminating  between  nuclear  explosions  and  earthquakes  and 
can  contribute  a  great  deal  to  understanding  regional  tectonics. 

Waveform  modeling  offers  the  best  hope  for  constraining  small-magnitude,  shallow  seismic  events  at  regional 
distances  with  sparse  observations.  Ideally  one  would  be  able  to  match,  and  thereby  determine  the  origin  of,  the 
time  and  amplitude  of  each  arriving  wave.  Such  a  match  is  an  unrealistic  goal,  in  general,  due  to  approximations 
required  for  tractability  in  modeling  algorithms,  to  an  incorrect  or  inadequately  precise  understanding  of  the  crust 
and  upper  mantle  in  most  regions,  and  to  poor  constraints  on  focal  mechanisms  and  source  time  functions.  While 
approaches  such  as  finite  difference  methods  offer  hope  for  accurate  and  computationally  feasible  three-dimensional 
(3D)  modeling  in  the  relatively  near  future,  to  be  useful  they  will  require  far  more  precise  and  accurate  velocity 
models  than  are  currently  available.  Assuming  that  accurate  3D,  or  even  2D,  modeling  is  currently  out  of  reach,  the 
question  that  arises  is  whether  ID  modeling  methods  can  be  used  in  a  strategy  that  minimizes  the  effects  of  laterally- 
varying  structure  and  focal  mechanism  to  constrain  the  focal  depths  of  seismic  events.  Presenting  and  evaluating 
such  a  strategy  are  the  goals  of  our  research.  Rather  than  matching  direct  and  reflected  pulse  shapes,  as  do 
Goldstein  and  Dodge  (1999)  for  larger-magnitude  (mb>=4.5)  and  more  distant  events,  we  seek  to  match  some  gross 
characteristics,  such  as  the  relative  arrival  times  of  a  series  of  waves,  as  well  as  possible. 

A  strategy  for  producing  ID  models  that  are  sufficiently  accurate  to  constrain  focal  depths  of  small  events 

Phases  that  arrive  near  the  direct  SV  phase,  including  Sp  (converted  at  the  base  of  the  Moho),  SsPmP,  and  shear- 
coupled  PL  (SPL)  waves,  collectively  sample  the  Earth’s  crust  and  upper  mantle  at  oblique  angles  and  therefore 
have  the  potential  to  produce  an  accurate  lateral  average  of  structural  properties  than  teleseismic  P  waves.  SPL 
waves  essentially  mimic  the  propagation  characteristics  of  regional  PL  phases,  with  the  important  difference  that  the 
number  of  events  available  for  modeling  is  often  greater  for  relatively  aseismic  regions,  since  sources  are  located  at 
teleseismic  distances.  SPL  waves  are  sensitive  to  crust  and  upper  mantle  structure,  including  seismic  velocity 
gradients,  Vp/Vs,  impedance  contrast  across  the  Moho,  and  layer  thicknesses. 

The  conventional  S  phase  (Figure  1)  is  the  initial,  relatively  sharp  and  pulse-like  arrival  that  signals  the  beginning  of 
a  wavetrain  with  generally  longer  periods  and  normal  dispersion.  The  particle  motion  associated  with  the  S  phase  is 
rectilinear  and  all  three  components  of  motion  are  in  phase.  The  dispersive  wavetrain,  that  follows  S,  exhibits 
prograde  elliptical  particle  motion  that  is  confined  to  the  vertical  plane.  Oliver  (1961)  named  this  wavetrain  “shear- 
coupled  PL’’  because  it  is  analogous  to  the  PL  wavetrain,  which  appears  between  P  and  S  arrivals  at  regional 
distances.  Oliver  (1961)  presented  a  theory,  based  on  the  observed  group  and  phase  velocity  of  SPL  that  explained 
the  phase  as  coupling  between  S  and  the  fundamental  leaking  mode  of  Rayleigh  waves  in  the  crustal  waveguide. 
According  to  this  theory,  shear  energy  generated  by  an  earthquake  (or  explosion)  travels  through  the  Earth’s  mantle 
as  a  body  wave,  whereupon  it  impinges  upon  the  Moho.  Afterward  a  portion  travels  through  the  crustal  waveguide 
as  trapped  P-waves  and  leaky  SV-waves  (Figure  2).  The  only  difference  between  a  PL  phase,  which  is  observed  at 
regional  distances  from  a  source,  and  SPL  phases,  which  are  observed  at  teleseismic  distances,  is  that  SPL  is 
generated  by  a  shear  wave  impinging  upon  the  Moho  at  regional  distances  from  the  observing  station.  In  addition  to 
producing  SPL  as  it  impinges  on  the  Moho,  a  portion  of  the  incident  S  wave  converts  to  P  as  well,  which  then 
travels  through  the  crust  to  arrive  at  the  station  as  a  precursor  to  S  (Figures  1  and  3).  This  phase  is  called  Sp,  and  it 
has  been  used  to  model  the  crust  by  Jordan  and  Frazer  (1975).  Its 
sampling  is  much  more  localized  to  the  station  than  is  SPL’s, 
making  its  sensitivity  less  representative  of  the  broader  region  and 
more  similar  to  that  of  P-coda  receiver  functions.  SsPmP  arrives  at 
the  base  of  the  crust  as  a  shear  wave,  travels  upward  through  the 
crust  as  a  shear  wave,  converts  to  a  P  wave  at  its  surface  reflection 
and  bounces  once  off  the  Moho  as  a  P  wave.  Langston  (1996), 
while  demonstrating  that  it  can  be  highly  useful  for  regional  crustal 
modeling,  showed  that  SsPmP  can  arrive  before  or  after  direct  S, 
with  either  larger  or  smaller  amplitude,  and  can  also  distort  the  S 
arrival  pulse.  We  will  attempt  to  simultaneously  model  S,  Sp,  and 
SsPmP  that  essentially  isolates  differences  to  the  P  structure  of  the 
crust,  for  data  collected  for  SPL  modeling.  Because  receiver 
function  methods  typically  deconvolve  the  vertical  seismogram, 
which  is  most  sensitive  to  compressional  wave  energy,  from  the 
radial  seismogram,  constraints  on  P  velocity  are  essentially 
sacrificed  in  order  to  obtain  clean  records  of  shear  phases.  The 
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Figure  1.  Typical  raypaths  for  S,  Sp, 
and  SsPmP.  These  body  wave 
phases  provide  constraints  on 
crustal  thickness,  P  and  S 
velocities  near  the  station. 


144 


26th  Seismic  Research  Review  -  Trends  in  Nuclear  Explosion  Monitoring 


data  we  propose  to  model  will  provide  a  valuable  check  of  receiver  functions,  in  that  they  constrain  the  bulk 
properties  of  the  crust— average  P  velocity  and  crustal  thickness  (Sp  and  SsPmP) — and  upper  mantle  (SPL).  If  the 
broad  model  search  turns  out  to  be  too  time-consuming  or  impractical,  a  fail-back  strategy  would  be  to  first  obtain  a 
basic  starting  model  by  modeling  Sp  and  SsPmP,  then  automating  the  SPL  modeling  using  the  Sp/SsPmP  model  as  a 
starting  point.  Since  reflectivity  is  a  full-waveform  method,  there  is  no  need  to  specify  which  phases  should  be 
included  (nor  need  we  identify  or  “pick”  specific  arrivals)  in  the  waveform  fitting  procedure.  While  the  waveform 
fitting  itself  will  be  automated,  a  great  deal  of  exploration  will  be  required  to  determine  optimal  window  lengths, 
frequency  content,  and,  perhaps,  variable  weighting  functions  for  different  portions  of  the  seismograms. 


vp* 


Figure  2.  Propagation  characteristics 
of  shear-coupled  PL  phases  (from 
Baag  and  Langston,  1985).  Note 
that  the  distance  of  propagation  of 
SPL,  and  therefore  its  sampling, 
depends  on  characteristics  of  the 
velocity  structure,  including  the 
slope  of  velocities  below  the  Moho 
and  attenuation.  The  earliest- 
arriving  and  largest-amplitude 
SPL  waves  are  those  that  have 
converted  from  S  nearest  to  the 
station,  so  our  modeling  will 
weight  local  sampling  more  highly 
than  distant  sampling  but  the 
wavefield  still  averages  structure 
laterally. 


Excitation  of  SPL-v/aves  with  distance  and  corresponding  period -phase  velocity  curve.  aN, 
N:  velocities  of  P -  and  £-waves  at  the  Moho,  T:  period,  V^\  phase  velocity 


Figure  3.  Depending  on  Earth  structure  and 
an  earthquake's  radiation  pattern,  the 
phases  Sp,  S,  SsPmP,  and  SPL  may 
appear  prominently  on  the  radial 
component  seismogram  at  distances 
between  30°  and  75°. 


Figure  4.  Statistics  of  an  example  run  to  compute  complete 
seismograms  for  frequencies  from  0  to  0.5  Hz  (number 
of  frequencies  =  2000,  number  of  layers=192,  number 
of  ray-parameters=l 000)  on  the  DEC  Alpha  cluster. 
Note  the  near  linear  speedup  of  the  algorithm  with  the 
increase  in  the  number  of  processors. 
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We  are  evaluating  the  usefulness  of  S,  Sp,  SsPmP  (Figures  1  and  3),  and  shear-coupled  PL  (SPL)  phases  (Figures  2 
and  3)  for  modeling  crustal  and  upper  mantle  structure  using  real  and  synthetic  data,  developing  a  waveform 
inversion  technique  based  on  a  novel  implementation  of  the  reflectivity  method  and  global  optimization  algorithms, 
and  applying  this  method  to  data  recorded  in  China.  We  have  made  substantial  progress  in  speeding  up  the  synthetic 
seismogram  computations  to  the  point  where  a  global  optimization  is  feasible.  The  reflectivity  calculation  involves 
computation  of  reflectivity  matrices  for  a  stack  of  layers  as  a  function  of  ray  parameter  (or  wavenumber)  and 
frequency.  The  computation  of  reflectivity  responses  for  different  ray  parameters  and  frequencies  are  completely 
independent  of  each  other.  We  took  advantage  of  this  independence  to  develop  a  reflectivity  code  that  runs  on 
parallel  computer  architectures.  Our  code  loops  over  ray  parameters,  i.e.,  to  each  node  we  assign  a  certain  number 
of  ray  parameters  to  compute.  At  the  end,  the  master  node  assembles  the  partial  responses  and  performs  the  inverse 
transformation  to  generate  synthetic  seismograms  at  the  required  azimuths  and  distances.  We  used  MPI  for  message 
passing  and  ran  our  code  on  a  PC  cluster  consisting  of  16  nodes;  each  node  is  a  660MHz  alpha  processor  with  8MB 
cache  and  1GB  of  RAM.  A  Myrinet  interconnect  is  used  to  communicate  between  nodes. 

Figure  3  shows  synthetic  seismograms  computed  using  our  parallelized  reflectivity  code  for  a  distance  of  50°  and  an 
earthquake  at  600  km  focal  depth.  Since  the  reflectivity  algorithm  is  “embarrassingly  parallel”  in  that  the  response 
for  each  frequency  or  ray  parameter  can  be  computed  on  independent  processors,  without  communication  between 
processors,  on  a  parallel  machine  computation  speed  increases  nearly  linearly  with  the  number  of  processors  (Figure 
4).  In  a  side-by-side  comparison  for  various  distances,  source  depths,  and  model  complexity,  our  code  matched 
results  of  the  Fuchs-Muller  code  very  well. 

Our  modeling  method  retains  the  time  and  cost  advantages  of  P-coda  receiver  function  methods  but  which  uses 
types  of  data  that  are  more  appropriate  for  nuclear  monitoring  purposes:  Shear-coupled  PL  phases  (SPL),  Sp  phases 
converted  at  the  Moho,  and  SsPmP.  SPL  samples  the  crust  and  upper  mantle  in  the  vicinity  of  a  station  most 
broadly  compared  to  Sp,  SsPmP  and  P  and  it  emulates  the  propagation  of  regional  phases,  which  reflect  at  more 
oblique  angles  (or  are  refracted  by  these  layers)  than  are  the  more  steeply  arriving  body  phases  typically  used  in 
receiver  function  modeling.  In  short,  because  of  the  data  they  use,  the  models  produced  by  receiver  function 
methods  may  be  inadequate  for  the  purposes  of  nuclear  monitoring.  These  latter  phases  sample  only  a  narrow  cone 
beneath  the  station  (e.g.,  Zhao  and  Frohlich,  1996).  Modeled  simultaneously  (where  they  exist),  SPL,  Sp,  and 
SsPmP  offer  the  potential  for  producing  azimuthally-dependent  structural  models. 

We  are  pursuing  this  strategy  because  efforts  to  determine  the  locations  of  small,  regional  seismic  events  are 
hampered,  in  most  parts  of  the  world,  by  insufficient  knowledge  of  the  crust  and  upper  mantle.  Also,  while  focal 
depths  are  often  a  highly  useful  discriminant  between  explosions  and  earthquakes,  their  determination  is  quite 
sensitive  to  crustal  structure.  The  most  encouraging  approaches  to  determining  focal  depths  require  precise 
modeling  of  seismic  waveforms,  particularly  for  small-magnitude  events,  in  which  travel  time  picks  are  relatively 
more  prone  to  errors  than  for  larger  events.  Yet,  a  precise  modeling  of  waveforms  at  regional  distances  requires  an 
accurate  model  of  the  crust  and  upper  mantle  along  the  propagation  path  between  source  and  receiver. 

RESEARCH  ACCOMPLISHED 

Forward  Modeling 

In  order  to  explore  the  sensitivity  of  each  phase  we  computed  synthetics  for  a  variety  of  distances  using  different 
crust  and  upper  mantle  models,  including  PREM  (Dziewonski  and  Anderson,  1981),  SNA  (Grand  and  Helmberger, 
1986),  TNA  (Grand  and  Helmberger,  1986),  ECH  (Zhao  et  al.,  1991),  and  WCH  (Beckers  et  al.,  1999).  Differences 
between  these  results  point  to  sensitivity  on  the  parts  of  these  phases  to  distinct  parts  of  the  models’  structure. 

Figure  5  shows  an  expanded  comparison  of  waveforms  produced  using  the  different  models  for  a  single  epicentral 
distance. 

Figure  6  shows  differences  in  waveform  features  to  be  expected  from  structural  variations  within  China  alone.  The 
models  at  upper  right  were  produced  by  Mangino  et  al.  (1999)  by  modeling  receiver  functions  for  structure  beneath 
stations  of  the  China  Digital  Seismographic  Network.  The  models  differ  primarily  in  the  lower  crust.  One  can  see 
from  the  synthetics  at  left  that  these  model  differences  are  reflected  by  differences  in  the  S  coda,  largely  in  SsPmP 
and  SPL,  for  the  distance  range  30°-50°.  These  synthetics  were  produced  for  a  600-km  deep  event.  This  suggests 
that  the  relative  timing  of  SsPmP-S  and  the  phase  and  amplitudes  of  SPL  phases  are  sensitive  to  the  lower  crust. 
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Figure  5.  Comparison  of  synthetic 
seismograms  computed  using 
various  Earth  models  for  the 
same  deep,  double-couple 
source  at  an  epicentral 
distance  of  30°.  Differences  in 
the  arrival  times  and 
amplitudes  of  S,  Sp,  SsPmP 
and  SPL  indicate  those 
phases'  sensitivity  to  structure 
near  the  receiver. 


Inverse  Modeling 

Our  modeling  process  is  controlled  by  a  global  optimization  algorithm  called  Very  Fast  Simulated  Annealing 
(VFSA)  (e.g.,  Sen  and  Stoffa,  1995).  Simulated  annealing  (SA)  is  analogous  to  the  natural  process  of  crystal 
annealing  when  a  liquid  gradually  cools  to  a  solid  state.  The  SA  technique  starts  with  an  initial  model  mo,  with 
associated  error  or  energy  E(mo).  It  draws  a  new  model  mnew  from  a  flat  distribution  of  models  within  the  predefined 
limits.  The  associated  energy  E(mnew)  is  then  computed  and  compared  against  E(mo).  If  the  energy  of  the  new  state 
is  less  than  the  energy  of  the  initial  state,  the  new  model  is  accepted  and  it  replaces  the  initial  model.  However,  if  the 
energy  of  the  new  state  is  higher  than  the  initial  state,  mnew  is  accepted  with  the  probability  of 
exp((^(m^  )-  £(m0 ))/  T),  where  T  is  a  control  parameter  called  temperature.  This  rule  of  probabilistic  acceptance 
(called  the  Metropolis  rule)  allows  SA  to  escape  local  minima.  The  process  of  model  generation  and  acceptance  is 
repeated  a  large  number  of  times  with  the  annealing  temperature  gradually  decreasing  according  to  a  predefined 
cooling  schedule.  A  variant  of  SA,  called  Very  Fast  Simulated  Annealing  (VFSA)  speeds  up  the  annealing  process 
by  drawing  new  models  from  a  temperature  dependent  Cauchy-like  distribution  centered  on  the  current  model.  This 
change  with  respect  to  SA  has  two  fundamental  effects.  First,  it  allows  for  larger  sampling  of  the  model  space  at  the 
early  stages  of  the  inversion  (when  “temperature”  is  high),  and  much  narrower  sampling  in  the  model  space  as  the 
inversion  converges  and  the  temperature  decreases,  while  still  allowing  the  search  to  escape  from  local  minima. 
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Second,  each  model  parameter  can 
have  its  own  cooling  schedule  and 
model-space  sampling  scheme. 
VFSA  therefore  allows  for 
individual  control  of  each 
parameter  and  the  incorporation  of 
a  priori  information.  The  model  is 
parameterized  in  terms  of  layers,  in 
which  Vp,  Vs,  density,  and  layer 
thickness  are  free  parameters. 
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Figure  6.  Synthetic  seismograms  computed  for 

various  receiver  function  models  of  the  crust  and 
upper  mantle  beneath  China  (models  from 
Mangino  et  al.,  1999).  The  primary  differences 
between  the  seismograms  are  in  the  S  coda, 
including  SsPmP  and  SPL,  which  indicates  these 
phases  are  highly  sensitive  to  distinguishing 
structural  details  of  these  models. 
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In  seismic  inversion,  more  than  one  model  can  often  explain  the  observed  data  equally  well  and  trade-offs  between 
different  model  parameters  are  common.  It  is  therefore  important  not  only  to  find  a  single,  best-fitting  solution  but 
also  to  find  the  uncertainty  and  level  of  uniqueness  of  that  solution.  A  convenient  way  to  address  these  issues  is  to 
cast  the  inverse  problem  in  a  Bayesian  framework  (e.g.,  Tarantola,  1987;  Sen  and  Stoffa,  1995)  in  which  the 
posterior  probability  density  function  (PPD)  is  the  answer  to  the  inverse  problem.  “Importance  sampling”  based  on  a 
Gibbs  sampler  or  a  Metropolis  rule  (Sen  and  Stoffa,  1998)  can  be  used  effectively  to  evaluate  the  necessary  multi¬ 
dimensional  integrals  and  to  estimate  PPD,  posterior  mean,  covariance  and  correlation  matrices.  The  posterior 
covariance  and  correlation  matrices  quantify  the  trade-off  between  different  model  parameters.  Sen  and  Stoffa 
(1995)  showed  that  multiple  VFSA  runs  with  different  random  starting  models  could  be  used  to  sample  models  from 
the  most  significant  parts  of  the  model  space.  This  “poor  man’s”  importance  sampling,  which  is  computationally 
efficient,  results  in  estimates  that  are  fairly  close  to  the  values  obtained  by  theoretically-correct  Gibbs  sampling. 

Example  of  an  Application  to  Data  from  the  China  Digital  Seismographic  Network 

Figure  7  shows  the  result  of  modeling  radial  and  vertical  records  from  station  BJT  of  the  9/28/1994  event  in 
Indonesian.  Subjective  choices  must  be  made  in  the  modeling  about  which  parts  of  the  seismogram  to  try  to  fit.  In 
regions  with  strong  lateral  variations  or  anisotropy,  it  may  be  that  no  single  ID  model  will  predict  all  the  phases 
adequately.  Requiring  the  modeling  to  try  to  fit  each  phase  similarly  will  result  in  a  model  that  does  not  predict  any 
single  phase  well.  Given  the  sampling  produced  by  the  various  phases,  we  tried  to  prioritize  fits  in  the  following 
order:  SPL,  SsPmP,  SV,  and  Sp.  Given  the  broad  sampling  and  lateral  averaging  performed  by  SPL  we  believe  that 
the  single  1 D  model  that  fits  SPL  best  will  be  the  most  useful  for  locating  earthquakes  regionally.  SsPmP  also  has  a 
broad  sample  and  constrains  Vp,  so  it’s  second  in  priority.  This  priority  order  reflects  a  starting  point — an  initial 
preference — and  will  be  evaluated  as  a  strategy. 

To  implement  this  strategy  we  successively  try  to  fit  broader  time  windows  of  the  data  waveforms.  8  15  shows  the 
results  of  600  iterations  in  which  only  the  window  between  815  s  and  860  s  is  evaluated  for  correlation  at  each 
iteration.  The  BJT  receiver  function  model  of  Mangino  et  al.  (1999)  is  shown  for  comparison,  although  it  is  not  a 
“reference”  or  “starting”  model  in  any  sense.  In  our  procedure  we  specify  the  search  limits  (also  shown  in  Figure  8) 
and  the  first  candidate  model  is  chosen  randomly.  Figure  9  shows  an  attempt  to  fit  a  broader  window  that  includes 
SV  and  Sp.  Differences  between  the  two  final  models  shown  in  Figures  7  and  8  are  again  concentrated  in  the  lower 
crust,  between  30  and  50  km  depth.  Figure  9  shows  the  result  of  a  600-iteration  run  to  model  records  from  BJT  of 
the  1 1/15/1994  Indonesia  event,  in  which  the  time  interval  fit  is  restricted  to  include  SsPmP  and  SPL.  Figure  10 
shows  the  results  of  a  600-iteration  run  to  fit  a  broader  window,  as  in  Figure  8.  Figure  1 1  shows  results  for  a  third 
event,  also  recorded  at  BJT.  Figure  12  compares  the  models  produced  for  the  region  containing  BJT  with  Mangino 
et  al.’s  (1999)  receiver  function  model,  in  the  context  of  the  imposed  search  limits. 

CONCLUSIONS  AND  RECOMMENDATIONS 

For  the  optimization  example  shown  in  Figure  8,  depth  zones  A,  B,  and  C  in  Figure  13  correspond  to  parameter 
correlations  indicated  by  the  same  labels  in  Figure  14.  Note  that  the  thin,  shallow  layers  in  A  are  characterized  by 
significant  off-diagonal  cross-correlations  (trade-offs)  with  other  parameters.  This  indicates  that  the  data  do  not 
constrain  these  layers  well.  Also,  the  symmetric  “variances”  about  the  mean  (actually  1  standard  deviation)  do  not 
contain  the  best-fitting  model  at  this  depth,  pointing  to  the  non-linearity  of  the  inverse  problem.  We  defined  the 
prior  to  be  Gaussian,  so  the  computed  (posterior)  variances  are  also  symmetric,  but  the  fact  that  the  best-fitting 
model  lies  outside  this  distribution  suggests  that  the  true  distribution  is  skewed.  Zone  B,  in  contrast,  contains  a  low- 
velocity  zone  in  Vs  that  appears  to  be  well-constrained  (Figure  13).  There  is  no  corresponding  Vp  low-velocity 
zone. 

Results  for  Zone  C  suggest  that  the  data  do  not  place  strong  constraints  on  the  model  parameters  at  these  depths,  on 
the  one  hand.  However,  since  the  best-fitting  model  is  pegged  against  the  search  limits  for  both  Vp  and  Vs,  we  will 
widen  the  search  limits  and  perform  a  greater  number  of  VFSA  iterations  to  explore  this  zone  further. 
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Figure  7.  Optimization 
results  after  600 
iterations  for  fitting 
synthetics  (black)  to  BJT 
vertical  and  radial 
records  (red)  of  the 
9/28/1994  event,  using  the 
time  interval  815-860  s. 
Best-fitting  model 
(black),  search  bounds 
(blue),  and  receiver 
function  model  by 
Mangino  et  al.  (1999) 
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Figure  10.  Same  as  in  Figure  9  except  that 
the  interval  800-860  s  was  used  in  the 
waveform  fitting. 
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Figure  8.  Same  as  in  Figure  7  except  that 
the  interval  800-860  s  was  used  in  the 
waveform  fitting. 
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Figure  9.  Optimization  results  after  600 

iterations  for  fitting  synthetics  (black)  to 
BJT  vertical  and  radial  records  (red)  of  the 
11/15/1994  event,  using  the  time  interval  815- 
860  s.  Best-fitting  model  (black),  search 
bounds  (blue),  and  receiver  function  model 
by  Mangino  et  al.  (1999)  (red). 
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Figure  11.  Optimization  results  after  400  iterations 
for  fitting  synthetics  (black)  to  BJT  vertical 
and  radial  records  (red)  of  the  5/24/1994 
event,  using  the  time  interval  835-870  s.  See 
Figure  9  for  other  details. 


Figure  12.  Comparison  of  models  for  the 
crust  and  upper  mantle  in  the  vicinity  of 
station  BJT  resulting  from  waveform  fits 
to  various  data  (see  legend),  Mangino  et 
al.’s  (1999)  receiver  function  model  (red), 
and  the  simulated  annealing  search  limits 
(blue). 
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Figure  13.  Mean  (cyan)  and 
best-fitting  (red)  models 
and  variance  of  the 
optimization  run  shown 
in  Figure  8.  Search 
limits  imposed  by  the 
operator  at  the  outset 
are  shown  in  green. 
Variances  are  only 
shown  for  Vs  (left)  and 
Vp  (right). 
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Figure  14.  Parameter  correlations 
computed  for  the  optimization 
run  shown  in  Figure  15.  There 
are  four  independent 
parameters:  Vp,  Vs,  layer 
thickness,  and  density,  repeated 
in  that  order  for  12  layers.  The 
total  of  48  parameters  is  shown 
here  in  color  (left)  and  with  a 
modified  grayscale,  in  which 
zero  correlation  is  white  and 
perfect  positive  and  negative 
correlations  are  both  black 
(right). 
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